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ENERGETICALLY STABLE DISCRETIZATIONS FOR CHARGE 
CARRIER TRANSPORT AND ELECTROKINETIC MODELS 

CHUN LIU*, MAXIMILIAN S. METTit, AND JINCHAO XU* 

Abstract. A finite element discretization using a method of lines approached is proposed 
for approximately solving the Poisson-Nernst-Planck (PNP) equations. This discretization scheme 
enforces positivity of the computed solutions, corresponding to particle density functions, and a 
discrete energy estimate is established that resembles the familiar energy law for the PNP system. 
This energy estimate is extended to finite element solutions to an electrokinetic model, which couples 
the PNP system with the Navier-Stokes equations. Numerical experiments are conducted to validate 
convergence of the computed solution and verify the discrete energy estimate. 
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1. Introduction and Background. Charge transport refers to any physical 
process where charged particles interact through an electric field and are driven by 
an electromagnetic force in some way. These systems have been observed throughout 
the history of science; they are, naturally, the foundation of electric engineering; and 
they are commonplace in everyday devices and physical systems, from mobile phones 
to solar-cell batteries, and weather to biology. The relevant quantities and their 
relationships are often modeled mathematically by the Maxwell equations, which were 
first published in [?]. To this day, phenomena relating to the transport and interaction 
of charged particles provides a broad variety of research topics in the physical sciences, 
mathematics, and engineering. 

In our study, we focus on electrostatic systems, where magnetic forces are absent. 
Such systems arise in biological settings, studying semiconductor devices, and elec¬ 
trokinetic systems, where charged particles interact with charged fluids, to name a few 
examples. Our model for charge transport is described by the Poisson-Nernst-Planck 
(PNP) system of differential equations. The electric field is defined by Gauss’ law in 
Maxwell’s equations, and the flux of the charged particles are driven by processes of 
diffusion and drift from the electrostatic force, which traces back to Nernst [?] and 
Planck [?]. This model is valid for systems where charge carriers can be accurately 
modeled as point charges. 

These equations serve as the basis for modeling many devices, such as batteries 
[?, ?], semiconductor devices [?, ?, ?, ?, ?, ?], fluidic micro/nano-channels and mixers 
[?, ?, ?, ?], and biological ion channels [?, ?, ?, ?], to name a few. As a system 
of coupled nonlinear partial differential equations, the PNP equations lead to a rich 
source of problems for pde analysis, where the system and its modifications are studied 
to improve understanding of the existence, uniqueness, and stability of a solution 
[?, ?, ?, ?]. 

Due to the wide variety of devices modeled by the PNP equations (or some mod¬ 
ification thereof), computer simulation for this system of differential equations is a 
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major application. This has led to a great deal of literature focusing on numerical 
solvers for the PNP systems [?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?]. Providing a compre¬ 
hensive numerical analysis would require an energy estimate to establish the stability 
of the discretization, some notion of convergence of the computed solution to the 
true solution, and well-posedness of the discrete problem. Prohl and Schmuck carried 
out such an analysis for the PNP system [?] and the PNP system coupled with the 
Navier-Stokes equation [?] using a numerical scheme that employs the method of lines, 
a finite elements discretization, and fixed point iteration. In this work, a novel finite 
element discretization is used that uses a logarithmic transformation of the charge 
carrier densities, which yields several favorable properties, such as automatic positiv¬ 
ity of the solution densities and energetic stability of the numerical solution for both 
the drift-diffusion model and the electrokinetic model. 

In §2, we define the PNP equations, introduce the energy law corresponding to the 
PNP system, propose our discretization and prove an energy estimates fully discrete 
solutions of the PNP system. In §3, we provide a similar analysis for the PNP system 
coupled with an incompressible fluid, where the divergence-free property of the fluid 
plays an essential role in establishing stability; consequently, this section also contains 
a discussion of a discontinuous Galerkin approximation for the Navier-Stokes system, 
which is known to preserve this divergence-free property. In §4, some numerical 
experiments are carried out to validate the convergence properties of the numerical 
solver and the discrete energy estimate. Some closing remarks are given in §5. 

2. The Poisson-Nernst-Planck equations and its discretization. The PNP 

system models the interaction of > 2 ionic species through an electrostatic field. 
We denote the ion density of the species by pi > 0 and the electrostatic potential 
by (j). Let If C for d = 1, 2, or 3, and T be a positive and finite real number. Then, 
the PNP system is described by the initial value problem: 

N 

(2.1) -V • (eV(()) = 

(2.2) = -V . t. 

(2.3) Ji — D'i^ Pi QiPiPi^(l)^ 

in n X (0, T], where 

Pi{x,0) = p^fi{x), for a: e n, i = 




and 


N 


-V 


eV(()(a::, 0)) = ec'^qiPi,oix), 


for X G V,. 


The electric permittivity, e = Ci-eo > 0, represents the the vacuum permittivity con¬ 
stant, eoj and the material dependent relative permittivity, e^, which may be discon¬ 
tinuous in general. The electric permittivity measures the strength of the long-range 
(nonlocal) interactions of charged particles. The term, Cc, represents the elementary 
charge constant. The ionic flux for the i**' ion species is denoted by Ji and is defined 
in (2.3) by a model proposed by Nernst [?] and Planck [?], where Di > 0 is the dif- 


fusivity, Qi is the valence number, and pi is the mobility of the T*' ion species. This 
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model is reasonable when the charge carriers are sufficiently small (with respect to 
the length scale of the domain) to be accurately modeled as point charges. 

We assume the Einstein relation holds, so that we may write 




where this relation implies that equilibrium distribution of the charge carriers should 
follow a Boltzmann distribution. Here, is the Boltzmann constant and T° is the 
temperature, which is considered to be fixed for the purposes of this paper. For 
simplicity, we choose our units of measurement such that Cc = 1 and kbT° = 1, so 
that fii — Di. 

2.1. Boundary conditions. The boundary conditions are a critical component 
of the PNP model and determine important qualitative behavior of the solution. 
A detailed account of stability and existence for steady-state continuous and finite 
element solutions has reported in [?, ?]. For the time-dependent case, existence and 
stability for the continuous case has been established [?, ?]; this work is concerned 
with establishing the stability of finite element solutions for the time dependent PNP 
equations, so homogeneous no-flux conditions are considered for each ion species, 

(2.4) Di{\7pi + qipi\7(j)) ■ n = Ji ■ n = 0, on dH. 

For the Poisson equation, write a disjoint partition of the boundary: dfl = Pd U Pat U 
Pfl with 


(j) = dV on Po, 

(2.5) eV(j)-n = S on Pa,, 

eV(j) ■ n + K(j) = C on Pfl, 


where 5V^ S and C are given functions. The Dirichlet boundary condition models 
an applied voltage, the Neumann condition models surface charges, and the Robin 
condition models a capacitor, where k > k > 0. Without loss of generality, it is 
assumed in the case where 5V is constant, so that 5V = 0. The capacitance is 
required to be positive on P^;, k > 0, though one may take Pi? = 0 if no capacitor 
is to be modeled. Any combination of Dirichlet, Neumann, and Robin boundary 
conditions can be applied to (j) for the purposes of this paper, though the case of pure 
Neumann boundary conditions requires the an additional constraint, which can be 
taken to be t) dx = 0 ioi 0 < t < T so that (j) is uniquely defined. 


2.2. Computational difficulties of the PNP system. The PNP equations 
present several difficulties when computing approximate solutions. Firstly, it is a 
strongly coupled system of nonlinear equations, which requires an iterative lineariza¬ 
tion procedure to resolve the nonlinearities, such as a Newton-Raphson method or 
fixed point iteration. While fixed point iteration serves as a helpfuf toof in the anal¬ 
ysis of the PNP system, it is difficult to establish its rate of convergence, which is 
critical in the practice of computing solutions. Secondly, the Nernst-Planck fluxes 
given in (2.31 are often convection-dominated, which leads to several analytical and 
numerical difficulties, such as the positivity of the ion concentrations, pi > 0, and 
the numerical stability of a given discretization. There are several ways to overcome 
these issues: the first is to introduce some sort of upwinding scheme, such as the 
Scharfetter-Gummel scheme and the box method [?, ?, ?], or the edge-averaged finite 
element (FAFF) method [?, ?]. Another option is to introduce a nonlinear change 
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of variables, such as the Slotboom variables [?, ?, ?] or the quasi-Fermi variables 
[?, ?, ?, ?], which symmetrize the Nernst-Planck flux. 

In this work, a novel change of variables converts the convection-dominated 
Nernst-Planck flux into a nonlinear-diffusion flux, similar to the quasi-Fermi variables. 
As a matter of fact, this change of variables is directly related to the quasi-Fermi vari¬ 
ables, though the quasi-Fermi variables introduce a nonlinear coupling between the 


equations in the time derivative term in (2.2). 


2.3. Energy of the PNP system. For PNP systems satisfying no-flux bound¬ 


ary conditions on the ion concentrations (2.4), it is known that the ion concentrations 
satisfy the conservation of mass. 


/ pi{x,t)dx= / pifi{x)dx, 
Jn Jn 


for 0 <t < T. 


Furthermore, in the presence of homogenous Dirichlet boundary conditions on (j>, the 
stability of the solution to the PNP system is known [?, ?] to be given by the energy 
law 


( 2 . 6 ) 


dt 


^p,(logp* - 1)-h I da;-f ^ ^ rfsj 

N 

X! 1^ (log Pi + dx, 



where the functional. 


. N . 

/ - 1) + I |V(/>|^ dx-h / ^ ds, 


is referred to as the energy and 


N 

E 

' i=i 


D^p,, |V (logp* -h q^(|))\^ dx > 0, 


as the rate of dissipation. The physical relevance of the no-flux boundary conditions 
on the ion concentration and the no-voltage boundary conditions on (j) stem from 
the notion that the PNP system is energetically closed; that is, there is no direct 
input or output of energy at the boundary. However, the case of applying a Dirichlet 
boundary condition to (j) is critically important in the analysis of many electrostatic 
devices, such as semiconductors, protein nano-channels, and electrokinetic devices. 
Accordingly, one can show that systems with a Dirichlet boundary condition on (f 
still satisfy a similar energy law. 

The energy associated with this system takes an unusual form compared to those 
typically encountered in finite element analysis due to the presence of the logarithm; 
nevertheless, this identity establishes the stability of the system as well as prescribes 
its rate of energy dissipation. Our choice of variables is motivated by the energy law 
(2.6), which specifies the regularity of the solution: take 
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which is the subpace of the usual Sobolev space, and for the ion concentrations, 


Pi € W = s p : I — y 


p(logp — 1) dx < oo and / |Vlogp| dx<oo 


L ' JO JQ ) 

leading implicitly to a positivity condition for the ions concentrations. A log-transformation 
of the ion concentrations yields a more familiar space 


p = logp G ir = n £“(1^), 


and, furthermore, guarantees positivity of the ion concentrations, since p = > 0. 

2.4. Log-density formulation and its energy. The standard £ 2 (f^) inner- 
product is used 


{u,v) 


uv dx, 


and inner-products on the boundary are given by 


i,v)r = / 
Jr, 


’ ds. 


and {u, v)n is similarly defined on Ttv- 

Using the log-density variables, the PNP equations are written in their weak form: 
find G W with r]i,t{t) G £^(n) and (j) G "Hp such that 


N 


(eV(/), Vv) -b {K(j),v)rn, - ^q^{e^%v) = {C,v)tr + (S',-u)r„, 

i=l 

-b (Die'^'W(r]i + qi(j)),'\/w) = 0, 

for z = 1,..., and all v G U, zc G W, and all times 0 < t <T, where 

{r]i{-,0),w) = for all w e W, 

N 

(eV(/)(-,0), Vv) -b 0),v)rn = '^qi{e'^'^''°\v) + {C,v)ra + (S',u)r„, for all z; G U. 




The energy law written in these new variables takes the form 
(2-7) + I + 

N 

= - / T Ae''’ |V {r]i + qi(j))\‘^ dx. 

i=i 

2.5. The discrete formulation. Let Tzi be a triangulation or tetrahedralization 
of the domain. For the usual space of piecewise linear polynomials, 

Wh = ^Wh G v\r G for all r G Thj C 'H^(n), 
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and denote the nodal interpolation operator, Ih ■ —)■ Wu- When Dirichlet 

boundary conditions are imposed on the electrostatic potential, define the spaces of 
continuous piecewise linear finite element functions 


VhXo = € Wh 

Vh,o = ^Vh G Wh 


=Ih{5V)], 
Vh\Tn = 0 } ■ 


When Robin boundary conditions are imposed, the lumped boundary inner-product. 


{u,v)R,h= / Th{uv)ds. 

Jtr 

is needed to for theoretical purposes to preserve monotonicity of the discrete Poisson 
equation. For the time discretization, define a partition of the time domain, 


0 — to ^ ^ * * * ^ 


where Atj = tj — tj-i- 

The finite element solution to the PNP equations is defined using the above finite 
element spaces and an implicit time discretization defined on the time partition: find 
%h ^ and e VhXo satisfying 

N 

( 2 . 8 ) (eS/(j)\i^\s/vh^ + (e'^'^Kvh^ = {C,Vh)R,h + {S,Vh)rr,, 

(2.9) ^ic/i) -h ^w/i), 

for i = 1,...,TV and all Vh G Vh^, Wh G Wh^ and j = 1,... ,m. The initial condition 
is given by 

(2.10) {r]l^l,Wh) = {r]t,o,Wh), ioiaWwhGWh, i = l,...,N, 

N 

(2.11) -f {C, Vh)R,h + {S, Vh)rN for all Vh G Vhfi. 

2.6. A discrete m 2 Lximum principle. The presence of a nonzero Dirichlet 
boundary condition imposes additional constraints on the hnite element mesh in order 
to maintain a discrete maximum principle for (j)h. Two approaches for satisfying 
a discrete maximum principle are summarized in the following lemma. The first 
approach constrains the interior angles of the mesh so that the discrete differential 
operator is monotone, the second approach requires quasi-uniformity and a sufficiently 
refined mesh. 

Consider an element, r G T^. The term facet is used below to denote an element 
edge when d = 2, and an element face when d = 3. Let E be an edge (1-dimensional 
sub-simplex) of r. The d — 2 dimensional simplex in r that is opposite to the edge, 

E, is denoted by (In two-dimensions, |fc^| = 1.) The angle, 0'f., is the angle 
between the facets containing edge E. The average value of e on element t is given 
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by {€)r = edx/\T\. In [?], it was shown that the off-diagonal entries of the stiffness 
matrix corresponding to the vertices on edge E are given by, 


loe = 


^(e)r|fcs|cot6»^ > 0, 


where the summation '^^^e taken to be the summation over all elements t G Th 
containing edge E. In the case where e is constant, this condition simply requires 
7ft to be a Delaunay mesh. Using this identity, a necessary and sufficient condition 
is given for the Poisson matrix to be monotone, implying that it has a nonnegative 
inverse and, consequently, a discrete maximum principle. 

Lemma 2.1. Suppose that one of the following assumptions hold: 

i. On each edge E G Th, it holds that uje = cot 6'^ > 0, 

ii. The permittivity, e, is constant and 7ft is quasi-uniform and sufficiently refined. 
Then, the finite element approximation, 4>h,D S Vft,r_D defined by 

(eV(^ft,D, Vuft) -I- {4>h,D,Vh)R,h = 0 for all vh G 14.0, 


satisfies a weak maximum principle 


max 10ft.£)(a;)I < Coo max \Ih{SV)\, 


for some Coo 4 1 that only depends on e, n, and the shape ofTl. Furthermore, under 
condition (i), the bounding constant is given. Coo = 1- 

The proof of case {i) follows from the monotonicity of the discrete differential 
operator and the proof for case {ii) can be referenced from [?]. In these works, there 
is no imposed Robin boundary condition, though the mass lumping discretization of 
this term preserves the discrete maximum principle. 


2.7. A discrete energy estimate. For autonomous boundary conditions, 6V, S, 
and C, the stability of the finite element solution analogous to (2.6) is verified. 
Theorem 2.2. 


Suppose 


G VFft 


and 4>h'^ 


G Vft.ro satisfy equations (2.8 ) — 


(2.11) for i = 1,... ,m and that one of the assumptions in Lemma 2.1 is satisfied. 
Then, the mass is conserved for each ion species, 


( 2 . 12 ) 




for i = 1,... ,7V, j = l,...,m. 


and the energy estimate is satisfied. 


N 




i<i<m 




1 


N 


+'^Atj / At 




7=1 


2=1 


dx-\--J Ih{K{(l)](’) )ds 


dx 




V0 


[FI + <ii&) 

dx+t f ds + C,. 

*7 r R 


where Ci depends on the number of ion species, their initial masses, the electric 
permittivity coefficient, and Coo- In the cases of no Dirichlet boundary conditions or 
a homogeneous Dirichlet boundary condition on 0ft, the constant Ci vanishes. 
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Proof. To prove the conservation of mass for the ion species, choose w/j = 1 G Wh 
in equation (2.9) to show 


Jn 


■^i,h — gPi.h dx = 


(j) 

e'i.h. — d'i.h 

At, 


, 1^ + (Ae''-^V(?7g + Vl) = 0, 


which yields (2.12). This argument expectedly fails when Dirichlet boundary condi¬ 
tions are imposed on since 1 ^ Wh- This is not an artifact of the discretization, 
however, and is also the case for the continuous system. 

f 7 ) f 7 ) 

For the energy estimate, set Wh = dl h + ^ which is a valid choice for 

( 7 ) 

the test function since (/)), G C Wh- This gives 


(j) r,G'-l) 

_ e’fi.h 

At; 




Atj 


■rf’ 


/ 

= - 

Jn 


ifil + tiJA) 


dx, 


which is summed over i = 1,..., A, to get 


(2.14) ^ 


N / JJ1 

e''-.'* - e''-.'* (,) 

-^dlh 


2=1 


Atd 




E 


% 


At, 


> Vh 


N 


(J) 


= -Y. 

.-1 7n 




dx. 


The first terms on the left are bounded by using the convexity of the function /(p) = 
p(logp — 1) for p > 0, which can be used to show 

(Pj - Pj-l) logPj > PjOogPj - 1) - Pj-i(l0gpj_i - 1). 

This follows from /'(p) = logp, f"{p) = Ijp > 0, and Taylor expansion. Applying 

(j ) ( J — 1) 

this bound with pj = and Pj-i = , one obtains for each i 


(2.15) 


U) (i-i) 

_ e^i.h 


At) 


Atj ’ 




- 1) - (e<h F _ 


0 - 1 ) 


,0-1) 


Pi.h > 


At, 


( 7 ) ( 7 ) 

To bound the remaining term on the left side of (2.14), decompose (jr^ = (PjIq + 

(t>h,D, where G Vhfi and 4>h,D G Vh^o satisfies the steady differential equation 
subject to the interpolated Dirichlet boundary condition: 


(2.16) ie'^'fh,D,'^Vh) + {K(j)h,D,Vh)R,h = 0, ’Ph.olrD =J'hiSy), 


for all Vh G 14 0 - Write 
(2.17) 


N / 

(,) 

> Vh 


E^i 

2=1 


At,- 


N 


= E'^* 


i=l 


„0) „0-i) 


At, 




N 


-E'l* 

2=1 




At, 




D 


and bound the first term on the right by subtracting consecutive time-steps of (2.8) 
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and taking vu = 


N 


2 = 1 


(i) 

li^h _ 


At, 


- 


= e 


-aE- 




Atj 

kO) Y7J.U) 


^hfi ^h,0 ,0)^ 

At, ’‘^^.0^ 

rf)\ 




At, 


i?,/l 


(2.18) 


> 




2 Atj 

lU) Aj) 




' R,h 


RM 


2Atj 


where the second equality follows from adding and subtracting the term 

At-^[{e\/(l)h,Dy(t)hl) + {l^(t>h,D,<Ph}))R,h] 

and the definition of (j)h, D, (|2.16[ ). 

Combining (2.14)— ( |2.15 ) and ( |2.17 ) —(2.18) gives the bound 


(2.19) 




(j) 


7 O- 1 ) 




At, 


N 


<-E / 




2 = 1 ' 


(-)S+<i.d’’) 


dx, 


where denotes the discrete energy functional, 


r N 




dx+l^ f Ih{K{(t>^^'')^)ds. 

^ JTn 


The first term in (2.19) yields a telescoping sum; a Gronwall argument leads to 

m N ^ 


max 

l<j<m 


’O') 




U) 


5]At, ^ / Ae''-:- 

,=i i=i 


(??s+ 




N 

- + E ■ 


2=1 


To complete the proof, the conservation of mass bounds 

N N 

(2.20) - e^^'’A4>h,D) < 2||g,e''c'> 


II II 

< 2N{ T^&j^ \\q,A'^.^\\^ )Uh,D\\c^, 


■ l<i<N 


c^n)t 


„(0) I 


where ||£ 1 (q) is directly proportional to the ionic mass, determined by the 

initial condition. The estimate, 

\\4>h,D\\c<^ < Coo max \Ih{5V)\, 

x^To 
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follows from Lemma 2.1 In the case where 5V = 0 or Ld = 0, it is clear that (j)h,D = 0 
so that this term vanishes altogether. □ 

To conclude this section, one important remark is in order. The inequality of 
this energy estimate is a consequence of only two aspects of the discretization; first, 
the time discretization satisfies ( |2.15| ) and (2.18) only with an inequality, whereas the 
semi-discrete solution (continuous in time) satisfies these bounds with equality. The 
only other inequality in the proof of Theorem |2.2| is used to bound non-homogeneous 
Dirichlet boundary conditions. As a matter of fact, in the semi-discrete case with 
homogeneous Dirichlet boundary conditions (or no Dirichlet boundary conditions), 
the finite element solution satisfies the energy estimate with equality. 


3. Electrokinetics. Electrokinetic systems combine effects of electrostatic sys¬ 
tems coupled with incompressible fluid flow. The model equations studied here cou¬ 
ple the PNP equations with the incompressible Navier-Stokes (NS) equations. This 
system of equations models electrokinetic phenomena such as electroosmosis, elec¬ 
trophoresis, streaming potentials, electrowetting, and many other phenomena where 
charged particles and fluids interact [?, ?, ?, ?]. Some analysis for this system in the 
continuous case is given in [?]. The equations governing the electrokinetic system seek 
a solution, rji,..., rjpf, (f), u, and p, such that 

N 

(3.1) _v.(eV0) = ^g,e''% 

(3.2) = V-(Ae’'-V(77, + ft0)-e''‘u), i = l,...,N, 

N 

(3.3) Pf{ut + (u ■ V)u) + Vp = V • (2/j,e(u)) — ^ V^, 

i=l 

(3.4) V-?I=0, 


on fl X (0,r], where e(-) denotes the symmetrized vector gradient, 

e(w) = (Vu)^). 

The initial conditions for this system are 


N 


r/^(x,0) = -V • (eV(/)(a;,0)) = 


i=l 

u(x,0) = uo(x), p(x,0) = po(x), for X € fl. 

Equations (O a nd (|3.2[ ) come direct y from the PNP model, where an additional 
coupling term in (3.2) models a kinetic force from the fluid flow described by the 
Navier-Stokes equations. Equations (3.3) and (3.4) are the usual NS equations for an 
incompressible fluid. In equation (3.3), the electrostatic force models 

the effects of the PNP on the fluid. 

The boundary conditions considered for the PNP variables remain the same as 
the previous section (2.4)-(2.5). The Navier-Stokes boundary conditions are assumed 
to be some combination of no-flux and no-slip boundary conditions. 


U ‘ Tl — 0, on Tno—flux ^ cAl, 

{fis{u)n) ■ f = 0, on Tno-flux, 

U — 0, on Tno— slip = \ hno—flux- 
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Due to the incompressibility condition on the fluid velocity (3.4), the solution satisfies 


e(u) = Vtt, which is commonly used to represent the viscosity term in the continuity 


equation (3.3). This identity does not hold, however, for general s G ["^^(^2)]“, so one 


must take care when using the divergence theorem to write the PNP-NS system in 
weak form; namely, for s = 0 on Pno-sUp, 

-(V • (2/u£(m)),s) = {2ne{u),e{s)). 


In the special case when m, s = 0 on 90, the right side reduces to (^Vm, Vs). 
The corresponding energy law for this system is given by 


(3.5) 


dt 



K 

2 


\(j>\^ ds 


Q. 


k(w)l {log Pi 

i=l 






dx. 


The terms in the energy law relating to the NS variables are critically hinged on 
specific mathematical structures of the NS system: in particular, the divergence- 
free property of the fluidic velocity plays a significant role in the cancelation of the 
cross-terms between the PNP and NS systems. As a result, the discrete solution 
must satisfy the divergence-free property on every subdomain of fl. This can be 
accomplished in several ways, using higher order elements or locally discontinuous 
Galerkin (DG) approximations [?, ?], for example. For many practical applications, 
solutions using higher order elements may be prohibitively expensive to compute; the 
discussion below primarily considers DG approximations for the NS variables. 

To define the weak solution to the NS equations, let 


Q = £2{^), 

S = {s G ['H^{n)f I u • n = 0 on Pno-flux, m = 0 on Pno-siip}, 
H(div;n) = {s G [C2{n)f I V • s G £ 2 (^ 1 )}, 

■Ho(div; fl) = {s G 'H(div; D) | s • n = 0 on 912}, 

■HD,o(div; 12) = {s G 'Ho(div; 12) | s = 0 on P„o_siip}, 

"HoCdiv®; 12) = {s G 'Ho(div; 12) | V • s = 0 in 12}. 

The weak solution of the NS equations is {u,p) G S x Q satisfying 

Dt{u\ u,s) -\- A{u, s) + B{s,p) = {f, s), for all s G S', 

B{u, q) = 0, for all q £ Q, 

where f G [£ 2 (^ 2 )]'^ and 

O f 

Dt{w\ u,s) = pf{ut,s) + • V)u,s), 

A{u,s) = {2pe{u),e{s)), 

B{u,q) = -(V • u,q). 
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The well-posedness of the weak formulation can be demonstrated using Babuska- 
Brezzi theory [?], where a Korn inequality must be established, as in the following 
lemma, which comes directly from [?]. 

Lemma 3.1. Let LI C = 2,3 be a polygonal or polyhedral domain. Then, 
there exists a positive constant Ck (depending on the domain through its diameter 
and shape) such that 

(3.6) 1^1 < C'ic||e(^||o, for all s € S. 


3.1. The discrete formulation. Recall that Th denotes the finite element mesh 
on LI and let St denote the set of interior element facets. The broken £2 and TL^ inner- 
products and norms are defined in the usual way 

ip,<i)n=^ip,ci)r, \\q\\o,n = iq,q)n, and = (Vs, Vs)^f, 

reTh 

for p,q € C2iLl) and s e H^^Th) = {w € C2{Ll) \ v\r € 'H^(r) for all r € Th}- 

Let w e T0-{Th), s G \}V-{Th)Y, ^ ^ \hL^{Th)Y^'^ denote a scalar, vector, and 
rank-two tensor field, respectively. These fields are "H^-regular within each element, 
though inter-element continuity is not assumed. Fix e £ Eh, where e = r_|_ H t_. 
Denote the outward unit normal vectors of r+ and t_ by n+ and n_, respectively; 
the averages across e on internal facets are defined by 

{u>} = ^(u;+ -bw-), {s} = ^(s+-bs_), and {ct} = ^(cr+ -b cr_), 

and given by their traces on the boundary facets; the jumps across internal facets are 
given by 


|w] = ri;_|_n+ -b W-fi-, 

[s] = s+ (g) n+ -b s_ (g) n_, 


[s] = s+ • n+ -b s_ • n_, 
[cr] = cr+n+ -b O-fi-, 


and lie] = wfi, [^ = s • n, |^ = s (g) n, and [a] = an on boundary facets, where the 
subscripts on the functions are equipped with their natural meanings of restriction 
to the element t+ or t_. An inner-product defined over the inter-element facets is 
defined 



w{s)v{s) ds. 


Using the facet average and jump notation, the following identities are readily verified 
by direct computation: for s • n = 0 on dLl, 



and 
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To preserve the local divergence-free property of the fluid velocity, nonconforming 
hnite elements are useful for assigning degrees of freedom aimed at preserving this 
property instead of conforming to the continuous spaces. We require the finite element 
space for the pressure Q/i C Q and Sh C ['H^{Th)]‘^ H "HD.oldiv), where Sh S, in 
general. While it is not necessary that St conforms to S, several constraints are 
imposed on the finite element pair Sh x Qh to ensure well-posedness of the discrete 
problem. First, it is required that 


(3.7) V • 5,, C Qh, 

and, second, that there exists for each qh G Qh a corresponding Sh G Sh such that 

(3.8) 's/-Sh=qh and ||s;j|jo < cp||gft||o, 


where cp > 0 is a Poincare constant that depends on Si in general, but not on qh- 
Requirements (3.71 and (3.81 together imply that V ■ Sh = Qh- The final requirement 
for well-posedness is the existence of a local interpolation operator. If,- : S\r —t Shir, 
for each element t G Th such that 


(3.9) 


|n^^l,r < C'ol^l.r and ||(/- n^)^|o,T < 


with hr = diam(r). This local interpolation property is extended over Tk, to yield an 
interpolation operator, flg,^ : S ^ Sh- 

Some well-studied hnite element pairs satisfying (3.7)-(3.9) are the Raviart- Thomas 
elements, Brezzi-Douglas-Marini elements, and the Brezzi-Douglas-Fortin-Marini el¬ 
ements all of degree fc > 1. Furthermore, as all of these elements are div-conforming, 
they have continuous normal components across inter-element facets, which, loosely 
speaking, “reduces” the discontinuity of the hnite element space, requiring simpler 
penalty functions in the discontinuous formulation. This additional continuity also 
plays a role in the cancellation of the PNP-NS cross terms and is commented upon 
in the proof of Theorem |3.2[ 


The discrete formulation of the NS equations given by: hnd iuh\ph^) G ShxQh 


such that 
(3.10) 


Dh,t{ur;F,^\sh) + Ah{uX\sh) + Bh{sh,pr) = ^\u 


(j)\ _ _PjLtoM ^),4) + (/(tj),s)i), for all 4 e S'/j 


(3.11) 


Bh{uh\qh) = 0 , 


for all qh G Qh 


for j = 1 ,..., m, where initial conditions are dehned by projection of the initial values, 
.-r(o) _ TT.. p ^^(0) 


= nSftUo and {pk ,qh) = (po, qh) for all qh G Qh- 


The forms used to dehne the discrete solution are given by 

Dh,t{wh;uh,Sh) =-^{uh,Sh) - ^{{wh ■y)sh,Uh) + / {wh ■ ftr)['^ ■ Sh) ds, 

2 

Ah{uh,Sh) = {2psiuh),e{sh))j-^ - {2p{eiuh)},lshj)sk - {2pluhj,{e{sh)})sk 

+ « X! ■ [41 ds, 

eeSh 

Bh{uh,qh) = —(V • Uh,qh)Th- 
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These forms are quite standard in the DG literature, though some terms and impor¬ 
tant properties remain to be specified. 

The discrete kinematic derivative term, Dh,t : S’;! x 5/1 x —)■ K, is defined using 

the upwind fiux, given by 


lim Uhix - Swh{x)). 
5-S.0+ ' ^ 


This definition yields coercivity, summarized by the standard identity [?]: 


(3.12) 


Dh^t{wh;uh,uh) 


1 

2 



• "ral I [wh] I ds, 


where n denotes either unit normal vector to the facet e. 

The bilinear forms, Ah and Bh, are a standard description for a DG discretization 
of Stokes’ equations and are motivated by the definition of the weak derivative followed 
by applying the divergence theorem element-wise. The parameter, a > 0, penalizes 
discontinuities of the solution across element interfaces and must be chosen to be 
sufficiently large. 

Since the finite element space, Sh, is div-conforming (3.7), equation (3.11) implies 
that V • Uft, = 0 on each element t G Tt- Another useful property inherited from (3.7) 
is that all Sh G Sh have continuous normal components across element edges; namely, 
letting ft and t denote the normal and tangent unit vectors, respectively, on each edge, 
e G 77,, gives 


Sh{x) = {sh ■ n)n + {sh ■ t)t = Sh{x) + sl{x), 
and |s7)l = 0 on each edge. As a result, it holds that 


idXfZ 


J[shj ■■ ads = Jlsll : ads for any a G [77^(D)]" 

Using this result, the coercivity of the kinematic derivative term, Dh^t-, reduces to 
Dh,t{wh\Uh,Uh) = I \wh ■ n\\luljf ds 


e^Sh 


and, for 

Ah{uh,Sh) = {2ne{uh),e{sh))j-^ - {2n{e{uh)},lSh\)Sh “ {‘^t^luhlAs{sh)})sH 



where only the tangent components along the element facets are penalized, and Bh 
becomes 


Bh{sh, Qh) = (V • Sh,qh) for all Sh G Sh and qh G Qh- 
The energy norm for the discontinuous fluid velocity is defined by 


INIIig — l^i.Th + where ^ ^11 Ptlllo,e- 
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For any of the three finite element examples mentioned above, one can establish the 
II ■ lluG'Stability of the bilinear form, meaning that there exists a positive constant, 
7 , such that 

(3.13) Ahish, Sh) > tII^IIdg: for all 4 e Sh (see [?]). 

This stability, together with an argument using fixed point iteration [?], is used to 
verify the existence of a discrete solution for the NS equations using this DG scheme. 

3.2. The discrete electrokinetic system. Employing the discretization of 
the PNP system given in the previous section and the discretization of the NS system 
above, the discrete solution to the electrokinetic system is defined by the finite element 
functions G Wh, G 14,n,, and ^P^h) ^ Sh x Qh satisfying 


(3.14) 


{eV(j)l^\Vvh) + {K4>h^ ,Vh)R,h 

II 

Vh) 





F{C^^\vh) 

R,h 4 

- {S^^\vh)rN, 

(3.15) 






A4 ^ ’ 

Wh) + ( 

^) , Vwh^ 

II 

■4 

,Wh) 


(3.16) 






Dh,t{u^h^ 

■,u)^\sh 

) +Ah{u\l\sh) A Bh{sh,p‘'h^) 

A4 ^ ^ 

) ^h) 

N 

f^h'^ ,Sh) 

i—1 

(3.17) 


Bh{'ifh^,qh) 

= 0, 




for all Wh G 114, Vh G 14,0) 4, € Sh, Qh S Qh, and j = 1,..., m. Initial conditions are 
prescribed by 


(3.18) 

{'nth'^h) = {m,0,Wh), 

for all Wh G Wh, 


N 

,\/vh) 4 {K(j)^°\vh)^j^ = '^q^{e'^'•°,Vh), 

i—1 


(3.19) 

+ {C, Vh)R,h + (4, 

for all Vh G 14,0, 

(3.20) 

40 ) TT - 
K =ns,^uo, 


(3.21) 

{Ph^^lh) = {po,qh), 

for all qh & Qh, 


where IIs,, : S ^ Sh satisfies ( |3.9[ ) locally. 

The stability of the discrete solution of the electrokinetic system is given by the 
following theorem. 


Theorem 3.2. Suppose -qfl G Wh,<t>h'’ 


VhXo^^h^ ^ Sh,Ph^ satisfy equations 
(3.14-) — (3.11), where Sh x Qh is one of the stable Stokes pairs as described above. 
Furthermore, suppose the mesh satisfies one of the assumptions of Lemma 2.1 Then, 
the mass is conserved for each ion species. 


JJ1 


{x,t) 


dx 



for 1 = l,...,iV, j = l,...,TO, 
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and the energy estimate is satisfied, 


(3.22) 


max 




N 

ai4’>r+i:e''S(,s-i)+i 


\7(j) 


U) 


1 


dx+-J Ih{K{4>h ) ) ds 




i=i 


II -(J)||2 

^hh Wdg 


i=i 


< ^ y |4°T + £ e"™ (4°?! “ + 2 / ds +Cl, 


where Ci depends on the number of ion species, their initial masses, the electric 
permittivity coefficient, and Coo- Iit- the cases of no Dirichlet boundary conditions or 
a homogeneous Dirichlet boundary condition on (fh, the constant Ci vanishes. 

Proof. The proof of Theorem |3.2| closely follows that of Theorem |2.2[ in addition 


to (3.13) for the NS variables. The only remaining terms are the cross terms between 


the PNP and NS systems, which cancel due to the strong divergence-free property of 
and the continuity of the normal components across element edges. 


The conservation of mass follows from choosing Wh = 1 € Wh in equation (3.15) 


as in Theorem |2.2[ Following the argument in the proof of Theorem 2.2 exactly gives 


(3.23) 




O) 


7O-1) 




Atj 


N 


<-Y. p- 


„(j) 

Pli,h 


1 = 1 ' 


ifd +f + v(,“ + ,. 0 “)) 


1 = 1 


N 


U) 


= -Y. 0.6”- 


dx+Y^ 


tip 


where the discrete energy is recalled as 


r N 




(fc) 


dx+^ [ Xh{K[(j)^4T) ds. 

^ JTr 


Let C,\^^ = Y4i=i Since strongly divergence free, has a contin¬ 

uous normE 
continuous. 


uous normal component across inter-element facets, • n = 0 on 9f2, and is 


N 

(3.24) Yi:' 




i=l 


(4-",vcf) 


= (v-4’.4>)r + E / p’' 


reTh 


I dr 


hf^>-fi^(Y)ds = 0. 
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Then, combining ( |3.23| ) and ( |3.24[ ) provides the bound 
(3.25) 






N 


At, 


U) 


< - V / 

-tn 




-O') 


For the NS terms, it follows from (3.12) 
(3.26) 


(3.27) 


)-|f( 


At, 


1 


+ b H /14" •^i|[(4" )]| 


(7") /I k 
and, choosing = Pj^ in (3.17), 

(3.28) 


> 


Bi.{4'’.p1’) = 0. 


Pf 


eeSh ' 
rd) I 


2At/ll llo 


70-1)1 


Setting Sh = ^ in (3.16) and employing the bounds (3.13), (3.27)-(3.28) gives 


(3.29) 


Pf 

2At, 


I-0)112 

14 llo-11“ 


O-OiiSn 


■7||»if>llLs-f;«(02vi,“,4’>). 


i=l 


Adding (3.25) and (3.29) gives 


(3.30) 


fii4'’iis 


+4* + ((hA'i.o) 


4"“ 


0 - 1 ) 112 
h Ho 


+ ^h + {Ch ^\4>h,D) 


< -7 M 


■U) I 


\DG 


Atj 

N 

-E 

i=l • 


D,e^ 


Jj) 


( 4 )+ 9 . 4 ’) 


dx. 


A discrete Grdnwall argument and invoking (2.20) gives the energy estimate. □ 


4. Numerical experiments. This section presents some numerical experiment 
that verify the viability, efficiency, and accuracy of computed solutions defined the 
proposed discretization in the above sections. According to the discretizations in §§2- 
3, a system of nonlinear elliptic equations must be solved at each time step. While 
there are many approaches to solving such a system, two commonly used techniques 
for resolving the nonlinear behavior are fixed point iteration (often referred to as 
Gummel iteration in the semiconductor literature) and Newton methods. For pur¬ 
poses of analysis, fixed point iteration is a very important tool, as convergence can be 
verified for more general problems; however, as a practical matter, it is often difficult 
to establish the rate of convergence to the nonlinear solution for this approach. This 
practical difficulty motivates the use of a quasi-Newton method for the experiments 
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presented below, where the relative residual approaches zero super-linearly. The na¬ 
ture of the model equations raise many issues concerning the numerical solver, such 
as resolving nonlinearities, upwinding schemes to preserve numerical stability, and 
describing the solver for the arising linear systems. Accordingly, some of the details 
of the numerical solver are deferred to an upcoming publication [?]. 

It is important to mention that the linearized equations resulting from a Newton- 
type approach lead to systems of linearized pdes that are potentially convection- 
dominated. This leads to potential algorithmic difficulties in preserving stability for 
the computed solution; so, some form of upwinding must be implemented to ensure 
accuracy. The well-studied edge-averaged finite element (EAFE) method is proven to 
provide stable numerical solutions that do not display spurious oscillatory behavior 
[?, ?]. A point of emphasis here is that the nonlinear solution is stable, as verified 
by the energy estimates above, thought the sequence of linearized equations are not 
necessarily stable. 

A solver was implemented in C-f-l- that leverages some existing functionality of 
the FEniCS 1.3.0 [?] software package for generating systems of linear algebraic equa¬ 
tions corresponding to an elliptic pde. Here, the elliptic pdes are the linearized pdes 
coming from Newton’s method, with an EAFE approximation to improve stability. 
Once the systems of algebraic equations are constructed, the Fast Auxilary Space 
Preconditioners (FASP) software package [?] is used to efficiently solve the resulting 
systems. 

The first experiment presented here is designed to establish the rate of convergence 
for the PNP discretization at steady state; since the discrete solution is defined using 
a method-of-lines approach, this experiment also verifies the convergence rate of the 
presented numerical scheme to the solution of the nonlinear elliptic equation at each 
time step. Several PNP systems are solved, where the permittivity coefficient, e, is 
tested for decreasing values. Testing the solver for small values of e is important 
in many practical problems concerning semiconductors and biological applications, 
where this coefficient is on the order of 10“^ to 10“® after the system has been non- 
dimensionalized. For this experiment, the equation 

-V • (eV</.) = +/o, 

=v(e''W(r?i+V<^)) + /i, 

=V-(e''^V(772-V(^))+/2, 


is solved on the domain H = [—1,1] x [—I, i] x [—J, j], where /o, /i, /2 are chosen so 
that the solution, for x = {xi,X 2 ,X 3 ) £ H, is 


r]i(x) = e 2 ("^1 1), ri2{x) = e * (“i+i)^ and (j){x) =— 


2 sinh(a:i) 


=-i 


As this experiment is designed to test the numerical convergence to the nonlinear 
solution at steady state, Dirichlet boundary conditions are imposed at the ends of 
the domain xi = ±1. The iteration count for convergence to the nonlinear solutions 
(deter min ed by reducing the relative residual by a factor of 10“^°) are reported in 
Table 4.1 along with the 'H} semi-norm of the error, given by 


semi—norm = -I- — ■ 


It is clear that the Newton iterates converge in a reasonable number of iterations 
(fewer than 10 in all cases), which is encouraging for small values of e. Additionally, 
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e 

20 X 10 X 10 

40 X 20 X 20 

60 X 30 X 30 

80 X 40 X 40 

1 

7 

2.65 X 10"^ 

6 

6.67 X 10-^ 

6 

2.97 X 10-4 

6 

1.67 X 10-4 

10-2 

6 

3.81 X 10-3 

6 

9.80 X 10-4 

5 

4.38 X 10-4 

5 

2.47 X 10-4 

0 

1 

ill- 

5 

7.03 X 10-3 

5 

2.37 X 10-3 

5 

1.20 X 10-3 

5 

7.18 X 10-4 

10 "® 

5 

7.25 X 10-3 

9 

2.61 X 10-3 

9 

1.43 X 10-3 

9 

9.34 X 10-4 


Table 4.1 

The count of Newton iterates to decrease the initial residual by a factor of 10“^® and the 'H^ 
semi-norm of the error. 



Fig. 4.1. The logarithm of the error measured in the 'H^ semi-norm, plotted against the loga¬ 
rithm of the element diameter. The lines depict the log of the error for various values of e, where 
the thick line is a reference for linear convergence. 


the convergence rate is to be linear for all values of e in Figure 141] with respect to the 
mesh size. 

A second experiment validates that the energy estimate is satisfied. While this 
property is certainly true for the theoretical finite element solution to the nonlinear 
problem, it is important to verify that the numerical solution, computed by inexact 
iterative methods, preserves this property. For this experiment, the domain is = 
[—1,1] X [—TO, to] X [—TO, to]. The time domain is (0, 0.03], and a uni form ti me-s tep of 
length At = For this problem, we solve the system defined by (2.1)—(2.3), with 
N = 2, fti = fj .2 = Di = D 2 = 1, qi = 1, q 2 = —1, and e = with no-rfux boundary 
conditions imposed on the Nersnt-Planck equations and mixed homogeneous Dirichlet 
and inhomogeneous Neumann boundary conditions on 4>- 


</' = 0 , 
eV0 • n = 1, 
eV4> ■ n = —1, 

eVcj) ■ n = 0, 


for a;i = ±1, 
for ail < 0 and Xz 
for ail < 0 and Xz 


for a;2 = ± 


1 

10 ’ 


1 

10 ’ 

1 

“10’ 


or ail > 0 and ^3 “ ~Yo’ 
or ail > 0 and ^3 ~ yf' 


These boundary conditions model surface charges, of alternating charge, lining oppo¬ 
site sides of a channel along the a;i direction and electrode contacts at the ends of the 
channel. The experiment demonstrates for each time step that the discrete energy 
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Fig. 4.2. The difference, log{—SE) — log(A), plotted over the time domain until convergence. 


estimate, 


cij) _ cO-1) 

= Ji - ^ - < _ 


+ e''^-'*|V(772.„ - dx = -A«, 


is satisfied until the dissipation is below the tolerance of the nonlinear solver, A^-^^ < 10’ 
To clearly illustrate that the energy estimate is satisfies, Figure [redisplays the quan¬ 
tity log(—— log(A('l)), which is positive when the energy estimate is satisfied. 


5. Concluding remarks. In this paper, the energetic stability is established for 
the finite element solutions to the PNP equations and an electrokinetic model with an 
incompressible fluid, with a minor extension to the case of inhomogeneous boundary 
conditions on the electrostatic potential. This extension imposes some additional 
constraints on the finite element mesh so that a weak discrete maximum principle 
can bound the energy introduced to the system by this inhomogeneous boundary 
condition. 

This energy estimate for the finite element solutions mimics the energy law of 
the continuous solution to these models, where a logarithmic transformation is a key 
ingredient to establishing the stability for the electrostatic terms and the divergence- 
free property of the discrete solution to the NS equations is essential to the stability 
of the fluidic variables, as well as the cancellation of cross terms between the two 
models. Recall that this divergence-free property of the finite element fluidic velocity 
is a result of using a DG formulation of the NS equations. 

A mathematical justification for the convergence is a matter of future work, 
though the experiments in this paper numerically demonstrate that convergence is 
obtained for various values of the permittivity coefficient in the Poisson equation. 
Furthermore, the computed solution using a quasi-Newton scheme is also shown to sat¬ 
isfy the energy law for the PNP system. The numerical solver for the PNP equations 
(with and without coupling to the NS equations) will be described in an upcoming 
publication [?]. 





